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We explore new ways of regulating defect behavior in systems of conservation laws. Contrary to 
usual regularization schemes (such as a vanishing viscosity limit), which attempt to control defects 
by making them smoother, our schemes result in defects which are more singular, and we thus refer 
to such schemes as “irregularizations”. In particular, we seek to produce delta shock defects which 
satisfy a condition of stationarity. We are motivated to pursue such exotic defects by a physical 
example arising from dislocation dynamics in materials physics, which we describe. 


A remarkable feature of systems of conservation laws 
is the tendency of smooth initial data to give way to 
non-smooth defects after only finite time. Examples in¬ 
clude shocks in hydrodynamics, cracks and dislocations 
in crystalline solids, and traffic jams in continuum models 
of traffic flow. Such defects are important in determining 
the qualitative behavior of a system, but they pose var¬ 
ious theoretical and practical challenges, stemming from 
the fact that the time evolution of defects cannot be de¬ 
termined solely from the evolution of the ambient con¬ 
tinuum. 

The physical origin of this difficulty can be understood 
as follows: Continuum equations of motion are derived 
from more fundamental, microscopic laws by discarding 
the “high-order” terms which are irrelevant at macro¬ 
scopic scales. Thus, systems with different microphysics 
can yield the same continuum equations, and conversely 
the continuum equations alone do not specify the mi- 
crophysics. Defects are inherently microscopic, and thus 
their behavior cannot be inferred from the continuum 
laws. 

The ambiguity of defect behavior in systems of con¬ 
servation laws is reflected in the mathematical theory 
by non-uniqueness of solutions containing defects. To 
obtain uniqueness, one must impose additional condi¬ 
tions upon the solution, presumably corresponding to the 
microphysics not accounted for by the continuum equa¬ 
tions. Common conditions of this sort include entropy 
eonditions and any of various types of regularization. In 
the former case, inequalities are imposed in analogy with 
the thermodynamical principle that entropy must always 
increase, while in the latter case the conservation laws 
themselves are modified by the addition of small terms 
that serve to smooth out defects. We will return in more 
detail to both of these cases in the next section. 

For the most common kinds of conservation laws (no¬ 
tably, those of hydrodynamics), the “proper” choice for 
solutions with defects turns out to be relatively tame, 
with defects consisting only of jump discontinuities and 
being described by regularization with a vanishingly 
small viscosity term. Some more exotic systems feature 


different classes of defects which require other types of 
regularization; for example, superfiuids are knowrP^ to 
exhibit dispersive shoek waves which are described by a 
very small dispersive term. Still another type of defect 
(which is a primary subject of this paper) is a delta shoek, 
which consists of a delta function-like spike, possibly rid¬ 
ing atop a regular (viscous) shock. This latter type of de¬ 
fect is much more poorly understood than the former two 
and differs in some important regards: Firstly, there is 
no particular regularization associated with delta shocks 
(for some systems, delta shocks may even be contained 
in the vanishing viscosity solution^. Secondly, there 
are few physically-realizable systems which are currently 
known to exhibit delta shocks. 



FIG. 1: Shocks in Burgers’ equation. (Color online) 
Various types of shocks evolving from smooth, 
sinusoidal initial conditions, depending on different 
regularizations of the Hopf equation: viscous shock 
(blue), dispersive shock (green), and delta shock (red). 

The present paper was motivated by the discovery of 
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delta shocks in a system arising from materials physics— 
continuum dislocation dynamics (CDD)—which to our 
knowledge provides the best physical interpretation of 
delta shocks currently available. These delta shocks have 
presented substantial challenges to physicists interested 
in simulating CDD numerically, which has led to a de¬ 
sire for new numerical methods tailored to better han¬ 
dle such exotic defects. We explore in this paper (in a 
mostly informal way) perspectives on delta shocks which 
may be of use to practitioners interested in simulating 
them numerically. (We do not, however, explicitly con¬ 
struct any numerical schemes here.) In particular, we 
describe a condition of stationarity which forces the for¬ 
mation of delta shocks, and we construct a new kind of 
“irregularization” which also produces delta-like defects 
(though we do not rigorously prove these to be true delta 
shocks). 

The outline of this paper is as follows: In the first sec¬ 
tion we provide some background on classical shocks. In 
the second section we describe delta shocks and discuss 
various ways in which they can arise in systems of con¬ 
servation laws, and we introduce the condition of station¬ 
arity and the irregularization mentioned above. Then in 
the third section we describe in more detail the system 
which initiated the preceding work—continuum disloca¬ 
tion dynamics. 


I. BACKGROUND 

In this section, we review basic facts about systems of 
conservation laws and standardize our notation. A sys¬ 
tem of conservation laws (or simply “conservation law”) 
is a partial differential equation (PDE) of the form 


+ V • F{u) = 0, (1) 

where u : Q x M+ ^ (for D C M’^), u = u(x^t)^ is 
the unknown, and F : MF is a given flux function. 

The conserved quantity u we refer to as “mass” in the 
following discussion. In the case where u is scalar valued 
(m = 1 ), we refer to eq. 0 as a scalar conservation law. 

Our working example of a conservation law will be the 
Hopf equation (a.k.a. the “inviscid Burgers’ equation”). 


referred to as a weak solution. Several distinct notions 
of weak solution exist, such as integral solutions, viscos¬ 
ity solutions, variational solutions, and various types of 
regularized solutions, to name a fev®. 

Not all forms of weak solution are applicable to all 
PDE’s. The types of weak solution which are most rele¬ 
vant for systems of conservation laws are regularized solu¬ 
tions and integral solutions. In the former case, the PDE 
is modified by the addition of a small term that serves 
to smooth out the defects; for the Hopf equation, such a 
modification looks like 


Ut+ { ] =V d^u, 


(3) 


where G M controls the strength of the regularization. 
To recover a solution to the original system of conser¬ 
vation laws, we take a limit as 77 ^ 0. Note that, in 
principle, the form of the regularizing term should cor¬ 
respond to the microphysics of a given system—so for 
example, systems with viscosity have a second-order reg¬ 
ularization, r]d‘^u^ whereas inviscid systems (like plasmas 
or superfluids) typically have dispersive (third-order) reg¬ 
ularization, gd^u. Below, in eqs. ( IdpT ), we will con¬ 
sider a different form of regularization altogether, which 
nonetheless bears a familial resemblance to eq. 0 . 

The main alternative kind of weak solution, the inte¬ 
gral solution, proceeds instead by integrating the original 
PDE to remove the problematic derivatives. Integral so¬ 
lutions are not of direct interest to us here, but we point 
out one tangential relationship to our current discussion: 
Integral solutions are generally not unique^ and to select 
a proper weak solution requires additional microphysi¬ 
cal information. This information is usually specified as 
a relationship between the values of the solution u on 
either side of a defect; when the relationship takes the 
form of an inequality it is called an entropy conditior^^^ 
and when it is an equality it is called a kinetic relatioiW. 
The stationarity condition we propose below is similar to 
a kinetic relation, but with the generalization that the 
speed of the defect can also figure into the equality. 


II. DELTA SHOCKS 


Ut + 



= UtF uux = 0 , 


( 2 ) 


which is a scalar conservation law in one spatial dimen¬ 
sion (n = 1). The Hopf equation is one of the simplest 
examples of a PDE which develops discontinuities in fi¬ 
nite tim^. Such discontinuities and other non-smooth 
features appearing in solutions to PDE we refer to col¬ 
lectively as defects. In the presence of discontinuities in 
the derivatives appearing in eqs. ([T][^ are no longer 
well defined, and so we must generalize our notion of 
“solution” to a PDE in order to make sense of the con¬ 
servation law at later times. Such a generalization is 


Eor our purposes, a delta shock (also written “^- 
shock”) may be defined informally as any defect which 
has a finite amount of mass concentrated at a point. We 
note in passing that there are also meaningful notions of 
(^'-shocks (and -shocks, more generally) correspond¬ 
ing to defects containing derivati ves o f delta functionJ^. 
Our numerical solutions (Section H C[ Eigs[^and[^ ap¬ 
pear to contain an upward and downward-pointing sin¬ 
gularity - the sum being the strength of the delta shock 
and the difference a 6' component. In this section we 
shall derive conditions for the evolution of the strength 
of the delta shock component. 

Delta shocks are known to arise naturally in a number 
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of systems of conservation laws—for example, the system 

) 

) (4) 


Ut + 

Vt + (wv)a 


considered by Tan, Zheng, & Zhan^, and the system 

Ut + -v)x = 0 

‘ ^ =0 (5) 


Vt+ \ -u -u 


studied by Keyfitz & Kranze]!^. In the former case, the 
delta shocks have been shown to be vanishing viscosity 
limits of the regularized system 


Uf {u 

Vt + {uv)x = 0 . 


(6) 


The mechanisms which cause delta shocks to form in 
systems of conservation laws are not fully understoocP, 
though in some cases the necessity of delta shocks is ob¬ 
vious. We highlight two such cases here: Firstly, if a 
conserved quantity v is a derivative of another variable 
14, and the latter exhibits regular (viscous) shocks, then 
clearly v must contain a delta shock. For example, if we 
take a derivative of the Hopf equation (§ and set v := Ux^ 
we find that 

Vf + {uv^x ~ Vjxt T {vjUx)x ~ ^x T ^2^ ^ ~ 

" ( 7 ) 

which is precisely the latter half of system Q. Thus 
(provided the initial data for system Q satisfies v{x^t = 
0) = u'{x,t = 0)) we see that this system can be inter¬ 
preted as the evolution of the conserved quantity of the 
Hopf equation, along with that of its derivative. 

This interpretation is not the only available, nor is it 
necessarily the best. Indeed, the closely related systenP 

Ut + {F{u))x = 0 

Vt + {g{u)v)x = 0, (8) 


shock balance the mass of the “hole” left as the shock 
moves away. In symbols, this reads 

F{ui) - F{Ur) =a{ui- Ur) , (9) 


where ui and Ur are the values of u to the immediate 
left and right (resp.) of the shock (see EvanJ^for a more 
formal derivation of this relation). 

If we allow the possibility of delta shocks storing mass 
on the defect, it is now the combination of the hole behind 
the shock and the mass m of the delta which must bal¬ 
ance the incoming flux. The modified Rankine-Hugoniot 
relation is thus 


F{ui) - F{ur) = cr{ui 


. dm 


( 10 ) 


Note that both eqs. ^ and ( p!Q| ) are vector identities in 
general, and may be interpreted as specifying the shock 
speed a and delta mass m (or the derivative thereof) in 
terms of the local environment, given by ui and Ur. In 
particular, this implies that eq. is only satisfiable in 
general when u is scalar, since vector u leads to an over¬ 
specification in § of the single unknown cr. Eq. ( p!Q| ) 
has no such problem, since m has the same number of 
components as 14, so that there is always one more un¬ 
known than there are equations in (10). This shows that, 
except in the presence of high degeneracy so that eq. 
is soluble for a, we expect multicomponent conservation 
laws to require delta shocks in order to conserve mass in 
the vicinity of a defect. 

This provides our preferred explanation of the delta 
shocks in system the classical Rankine-Hugoniot 
lation ^ for this system reads: 


re- 


— ul = a {ui — Ur) 
UiVi — UrVr = CF {vi — Vr) 


a = 


a = ui F Ur 
UiVi — UrVr 
Vi - Vr 


, (11) 


which does not hold unless UiVr = UrVi (i.e. the necessary 
degeneracy for © to be satisfiable does not generally 
exist). A similar statement shows why the more general 
system also displays delta shocks. 


where F is an arbitrary smooth & convex function and 
g is an arbitrary smooth & increasing function, also dis¬ 
plays delta shocks, though they cannot in this case be in¬ 
terpreted as derivatives of regular shocks. A more general 
interpretation for the appearance of delta shocks in such 
conservation laws is as a consequence of incompatibility 
of defect motion and the flow of the conserved quantity^ 
which we explore in the next two subsections. 

A. Generalized Rankine-Hugoniot relation 

The Rankine-Hugoniot relation is a statement of con¬ 
servation of mass across a shock. Eor the case of a regular 
shock in one spatial dimension which moves with speed 
cr, mass conservation requires that the net flux into the 


B. Stationarity 

The Rankine-Hugoniot relations illustrate how multi- 
component systems can develop delta shocks when a de¬ 
fect’s motion cannot be simultaneously compatible with 
the flow of all components of the conserved quantity. We 
can also produce a similar effect even in scalar conser¬ 
vation laws by insisting that a defect move according to 
specified kinematics which are not compatible with the 
flow of the scalar conserved quantity. The simplest such 
example is to demand that a defect be stationary. This 
might be an appropriate kinematic constraint in, for ex¬ 
ample, an electrical network containing fuses, where a 
blown fuse is modeled as a stationary defect. A station¬ 
arity condition could also be physically appropriate for 
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the continuum dislocation dynamics example we will de¬ 
scribe below. This stationarity condition is analogous to 
a kinetic relatioiP, with the generalization that the speed 
of the defect can now figure in as well (rather than just 
the value of the local conserved quantity). 

As an illustration, we consider the Riemann problem 


u{x^ t = 0) 


{o 


X < 0 
X > 0 


( 12 ) 


for the Hopf equation (§. Mass flows into the defect at 
- = 0 with rate F(O-) - F(0+) = ^ • | • 0^ = 


and thus by the modified Rankine-Hugoniot relation (10) 
with cr = 0 we find that the mass of the delta at x = 0 
grows at a rate m = |. Everywhere else the value of u is 
uniquely specified by projection of characteristics. Hence 
a reasonable candidate solution for all times t > 0 is 


u{x, t) = u{x, t = 0) + ^ 5{x) (13) 

We note that this candidate solution is not an inte¬ 
gral solution in the usual sense, as the square of a delta 
function (needed to make sense of an integral solution 
to eq. ([^) is undefined by classical distribution theory. 
(However, it is possible to overcome this technicality us¬ 
ing generalized distribution^.) 


C. Regularizations for stationarity 

Our ad hoc construction of a solution ( p^ to the Hopf 
equation which satisfies the condition of stationarity does 
not lend itself well to numerical simulation. We thus find 
it advantageous to try to construct a “irregularization” 
of the Hopf equation which yields such a solution. Our 
attempted irregularization looks like: 

Ut+ =0, (14) 

where the smooth & bounded function f{e^Ux) satisfies: 

1. /(e,0) = l 

2. /(e,-oo) = 0 

3. /(0,-) = l. 

For example. 


Ut + 


\v? 


eui 


= 0 


X 


(15) 


above implies that / ~ 1, and the original Hopf equation 
is recovered (approximately). 

(This prescription for the irregularization is not as gen¬ 
eral as possible—for example, one also could add to the 
denominator of eq. (15) a non-vanishing polynomial in 
which would reflect some sort of self-interaction effect. 
However, the “minimal” irregularization of eq. (15) is suf¬ 
ficient to produce the defects we are interested in. More¬ 
over, the effect is not peculiar to this particular form: We 
have also experimented with the irregularization 



and seen the same qualitative behavior. ) 

For numerical work, it is sometimes helpful to also add 
a small regularizing term to ensure solutions stay rela¬ 
tively smooth—e.g. 


+ 




= n d^u, 


where k is typically 2 or 4. 


(16) 



FIG. 2: Irregularized Hopf equation. (Color online) 
Simulation of the irregularized Hopf equation ( [l^ 
starting with sinusoidal initial conditions, for two values 
of e at time t = 1 (simulated with a conservative upwind 
differencing). Note that the singularity has both a 
positive and negative peak. Only the sum contributes 
to the delta shock (as a weak limit); the difference can 
be described as a derivative of a (5-function. 


is such an irregularization. The intuition behind this is 
that when a shock develops in the Hopf equation, Ux 
diverges to —oc, sending the modified flux f{e^Ux) 
to zero. The quenching of the flux will hopefully cause 
a pileup of mass behind the defect. On the other hand, 
for small values of e and away from defects, condition 3 


Does this succeed in making the defect stationary? 
We b^in to answer this by examining numerical results. 
Fig. (0 shows simulations of eq. (16) for e = 0 & e = 10“^ 
on the unit interval with periodic boundary conditions, 
for parameters t = l,/c = 2,?7 = 10“^, and a grid spacing 
dx of .002. 
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Clearly the irregularization has substantially slowed 
the movement of the shock and thereby induced a delta- 
like spike to conserve mass. The extent to which this 
is a bona fide delta shock is a delicate question (we will 
thus cautiously refer to this feature as a delta spike). We 
mention three relevant considerations: 


1. To begin with, the effect is inherently non-smooth, 
as for smooth solutions of eq. (15), we can prove 
the following maximum principle: 


Proposition 1. A solution u to eq. on a eompaet 
domain satisfies 


max{u{'fi)}=max{u{'^0)} (17) 


for all times t > 0 while the solution remains smooth. 


We motivate this as follows: If we look at a curve 
y(t) such that attains a maximum at y{t)^ 

then we find 


dtu{y{t),t) 


— T 

= u^y- 

= - 

= 0 


\u^ 


1 + 


UxU 

-\- eui 


eUx ^xx'^ 


(1 


(18) 


Where the last expression vanishes because Ux = 0 
at the maximum y{t). This shows that the maxi¬ 
mum value of u is unchanging in time. (This ar¬ 
gument can be made more rigorous, avoiding the 
assumption that the maximum can be parameter¬ 
ized by a differentiable curve ^(t), but we will not 
do so here.) Note that a similar result holds also 
in the case of eq. (16) with a 2nd order (viscous) 
regularization, with the modification that the max¬ 
imum of u is then decreasing in time. 


2. Examining the delta spike closely reveals an inter¬ 
nal structure. To provide some perspective, we note 
first that eq. ( p!^ admits two relevant exact solu¬ 
tions, one of the form a{t) -1- b{t)x (i.e. a straight 
line, moving in time), and the other a stationary 
solution of the form 


a cosh 



(19) 


Experimentally, we find that the shape of the spike 
is well approximated by a profile of the form (19), 
as shown in fig. More precisely, the form of 
the delta spike consists of two partial coshes, one 
positive and one negative, separated by an abrupt 
discontinuity. Away from the defect, the solution to 
eq. (15) approaches the aforementioned linear solu¬ 
tion asymptotically in time. This eventual shape is 


analogous to the asymptotic “N-wave” profile of so¬ 
lutions to the regular Hopf equation. As discussed 
in the appendix, it is possible to derive analytic 
expressions for the parameters a(t), 6(t), a, (3 de¬ 
scribing this asymptotic form, and we find good 
agreement between simulation and these analytic 
values. 



EIG. 3: Cosh profile of delta spike. Comparison of a 
delta spike with a cosh profile of form (19). Note that 
the cosh profile shown here is not a best fit, but rather 
an analytically derived approximation, as discussed in 
the appendix. The delta spike comes from a simulation 
with parameters e = 2.5e — 6, dx = 2.5e — 4, and initial 
data u{x^A) = sin(27rx) -h .1. Again, the positive and 
negative spikes must be added together to give the net 
weight of the delta shock. 


An important feature of the cosh profile (19) is that 


its width scales like yA. To conserve mass, the 
height of the shock must therefore scale like 1 /yT. 
This suggests that, although the delta spike is not 
truly a delta shook for non-zero values of e (indeed, 
the cosh profile (19) certainly does not qualify as a 


delta shock), as e ^ 0 the spike does indeed con¬ 
verge to a delta function. Experimentally, this is 
what we seem to find (conditional upon some fi¬ 
nessing of the numerical methods), as discussed in 
the next point. 

3. Establishing numerical convergence of our simula¬ 
tions is tricky for two reasons: Eirstly, since we are 
looking for convergence to a delta shock, we can¬ 
not hope to have convergence in e.g. the L 2 sense. 
Secondly, the equation (15) is highly unstable (not 
surprisingly, since the irregularization is designed 
specifically to produce a delta shock, which may be 
viewed as a sort of instability). After trying various 
finite-differencing schemes, we found that the best 
results came from a type of conservative upwind 
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differencing, wherein u is differenced upwind but 
Ux is not. Explicitly, the semi-discrete formulation 
of eq. (14) is 


dtUi + — 0 


/ (e’ + ^i+l) > 0 

1 «|+i/ (e, {Ui + Ui+i) < 0 


( 20 ) 


where A is the discrete forward difference opera¬ 
tor {Aai = — ai). All figures presented herein 

were produced using this discretization. (Remark: 
for the regularized eq. ( p^ we performed an opera¬ 
tor splitting, evaluating the regularization term in 
Fourier space for stability.) Despite being the most 
successful numerical method we found, the stabil¬ 
ity properties of this discretization are not entirely 
satisfactory. In particular, we found that when 
the grid size became substantially smaller than the 
width of the delta spike (which is on the order of 
V^), the delta spike would become unstable and 
break apart into multiple, smaller delta spikes. To 
overcome this difficulty, we let the mesh size dx and 
the width of the delta yT go to zero together, keep¬ 
ing dxj y/t fixed (recall that it was our intention to 
let e ^ 0 anyhow, so as to recover a solution to the 
original Hopf equation). This technique appears to 
work, as seen in fig. The delta spike gets nar¬ 
rower and higher, while remaining in essentially the 
same place, as dx, e ^ 0 in this manner. 



0.525 0.530 0.535 0.540 


FIG. 4: (Color online) Close-up of the delta defects as 
dx, e ^ 0, r\j dx. 


Although our simulations seem to establish conver¬ 
gence to a delta shock, it is a very peculiar convergence 
wherein the strength e of the irregularization must van¬ 
ish together with the discretization length dx. The fact 
that the width of the delta spike must be kept close to 


dx suggests that there may be some sort of resonance be¬ 
tween the irregularization of eq. and the grid which 
drives the formation of the spike (especially in light of the 
maximum principle mentioned previously). Nonetheless, 
we can still hope to produce a meaningful weak solution 
for the Hopf equation in this way. We will not establish 
anything more rigorous at this time, though. 


III. PHYSICAL MOTIVATION: DELTA SHOCKS 
AND DISLOCATION WALLS 

Our interest in delta shocks is motivated by a tangible 
phys ical q uestion. Why do dislocations in crystals form 
wall j^^ l ^^ -^ Briefly, a crystal is a regular array of atoms. 
A dislocation line is a flaw in that array, such as the edge 
dislocation formed by the boundary of an extra plane 
of atoms, or the screw dislocation formiM the central 
line where planes form a ‘spiral staircase^. These dis¬ 
locations move to mediate plastic deformation when the 
crystal is bent, and form tangles that organize into wall¬ 
like structures called cell wall^^. In a continuum theory 
of dislocation dynamics, such walls must be described as 
delta shocks: their density scales as the inverse square 
of the lattice constant, which vanishes in the continuum 
limit. 

One c lass of continuum dislocation dynamics theo- 
riei ji5 | 20 | form such delta shocks. In these theories, 
the dislocations move with a common velocity J^. If one 
describes an incipient wall in the yz plane with a dislo¬ 
cation density that depends only on x, the dependence 
of this velocity on the dislocation density simplifleJi^ to 
the Hopf equation: 


=0 


( 21 ) 


When the Hopf equation forms a step-like shock, the en¬ 
trained dislocations pile up into a delta shock, forming a 
wall. 

The central physical question is how that wall should 
evolve after it forms. If we regularize the Hopf equa¬ 
tion (21) into Burger’s equation (i.e. viscous regulariza¬ 
tion) , the dislocation wall moves along with the Rankine- 
Hugoniot velocit y. This produces a self-consistent, sen¬ 
sible modeP^^, but one that is unsatisfying in two re¬ 
gards. First, the velocity of dislocation walls in crystals is 
determined by the microstructure of the walls and not by 
continuum properties. Secondly, generalizations of this 
continuum theory have differing velocities for the differ¬ 
ent components of the dislocation density - smearing the 
resulting walls. 

On the other hand, physically at some junctions where 
two dislocations intersect, the point of intersection can be 
pinned in place (so-called sessile dislocation junction^^. 
Dislocation walls can also be pinned by impurities that 
segregate to the boundar y. Th is provides the motivation 
for our study in section (HB) of the Hopf equation (21) 
with stationary shocks. 
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We should note that the behavior of both our model 
in two and three dimensions and experimental disloca¬ 
tion systems under stress is more complex than simple 
wall singularities. There they produce complex cellular 
structures, which in our models and some experiments 
form self-similar, fractal morphologies. Our theories in 
higher dimensions show clear analogieJ^ to the behav¬ 
ior of the Euler equation (the inviscid limit of Navier- 
Stokes). Nonetheless, these structures are complex, ram¬ 
ified wall-like entities, whose dynamics should be con¬ 
trolled by physics on the microscopic atomic scale, not 
by the continuum laws. 


IV. CONCLUSIONS 


We have seen that delta shocks arise naturally in a va¬ 
riety of systems of conservation laws, and their presence 
can often be traced to an incompatibility of the flow of 
the conserved quantity and the motion of defects. Al¬ 
though most of the existing examples of delta shocks do 
not have direct physical interpretation, the case of cell 
walls in crystals does appear to furnish such an example. 

A recurring theme in our presentation has been the im¬ 
portance of developing numerical techniques capable of 
handling defects with more exotic behavior than tradi¬ 
tional, viscous shocks. Ideally, one would like to be able 
to specify whatever defect behavior is believed appro¬ 
priate for a given system and have numerical techniques 
which respect that behavior. In lieu of such a very gen¬ 
eral approach to simulating conservation laws, we must 
instead focus on specific classes of defect behavior which 
we believe are important. To this end, we propose further 
study of the stationarity condition introduced above. 

We have concocted an “irregularization” ( IdpT ) of the 
Hopf equation which appears numerically to exhibit delta 
shocks, suggesting a possible avenue to development of 
the aforementioned numerical methods for delta shocks. 
Though we are far from showing that these apparent 
deltas are bona fide, we nonetheless think that the ir- 
regularization and the defects it yields are interesting in 
their own right and merit continued study. 


V. APPENDIX 


We present here some further numerical and analytical 
observations (mostly without details or proofs) for the 
irregularization 


Ut + 


1^,2 


■ e ut 


= 0 . 


X 


( 22 ) 


In particular, we sketch a systematic way to approximate 
the asymptotic form of solutions to this equation. 

Three useful characteristics of the regular Hopf equa¬ 
tion can be shown to also hold for smooth solutions of 
this more general equation: 

1. The maximum principle mentioned above. 


2. Zero-crossings of a solution to eq. (22) remain fixed 
in time. 


3. The total mass between any two zero-crossings is 
conserved. 


The second and third items above allows us to predict 
(analytically, approximately) the asymptotic form of so¬ 
lutions to eq. ( [ 2 ^ . We sketch this analysis as follows: 
Shocks tend to form on downwards-going zero crossings, 
and away from shocks the solution becomes linear at large 
times. The overall profile of the solution thus is linear ev¬ 
erywhere except at the shocks, where it is approximated 
by the cosh profile mentioned above. The stationarity of 
upwards-going zero crossings allows for analytical evalu¬ 
ation of the slope there, which asymptotically must equal 
the slope of the entire linear region of the solution. Ex¬ 
plicitly, the slope at a zero crossing Xzc satisfies the 
ODE 


0 = dtU^{xzc,t) + - -- - -^ 

2 \l^eu^{xzc,ty 






1 + eUx{xzc,ty 


(23) 


which is obtained by differentiating eq. ( |22| ) w.r.t. x and 
noting that all terms with a factor of u vanish. This ODE 
(23) can be solved analytically, providing the asymptotic 
form of the linear part of the solution. 

The third item above allows us to evaluate the mass of 
the delta, which is just the difference between the total 
initial mass between two zero-crossings and the eventual 
mass of the linear solution on the same region. Examin¬ 
ing the cosh profile of the delta mentioned above. 



we see that a must equal the height at which the delta 
begins, which can be determined from the linear solution. 
Thus a is determined, and then (3 can be determined from 
the known mass of the delta. 

Stringing together the above arguments allows a piece- 
wise description of the solution across the entire domain. 
This analytic, approximate solution is found to agree rea¬ 
sonably well with simulation—for instance, the cosh pro¬ 
file shown in fig. ^ uses parameters calculated in this 
way. 


VI. ACKNOWLEDGEMENT 

We thank Randall LeVeque for helpful comments and 
guidance on hyperbolic conservation laws and numeri¬ 
cal methods. We are especially indebted to Alexander 
Vladimirsky for his exceptionally insightful suggestions 
for the production of this paper. This work was sup¬ 
ported by US DOE/BES grant DE-EG02-07ER46393. 









^ P. Lax, Hyperbolic Systems of Conservation laws and the 
Mathematical Theory of Shock Waves, Society for Indus¬ 
trial and Applied Mathematics, 1973 
^ P. Lax, Hyperbolic Systems of Conservation Laws //, Com¬ 
munications on Pure and Applied Mathematics, 1957 
^ P. G. Le Floch, Hyperbolic Systems of Conservation Laws: 
The Theory of Classical and Nonclassical Shock Waves, 
Birkhuser Verlag AG, c/o Springer GmbH & Co, 2002 
^ L. Evans, Partial Differential Equations, second ed., Amer¬ 
ican Mathematical Society, 1998 
^ D. C. Tan, T. Zhang, T Chang, and Y. X. Zheng, Delta- 
Shock Waves as Limits of Vanishing Viscosity for Hyper¬ 
bolic Systems of Conservation Laws, Journal of Differential 
Equations, 1994 

® J. E. Colombeau, New Ceneralized Functions and Multipli¬ 
cation of Distributions , Elsevier Science Publishers B.V., 
1984 

^ B. L. Keyfitz and H. C. Kranzer A Viscosity Approxima¬ 
tion to a System of Conservation Laws with No Classical 
Riemann Solution ’’Nonlinear Hyperbolic Problems,” Lec¬ 
ture Notes in Mathematics, Springer-Verlag, 1989 
® B. L. Keyfitz, Conservation Laws, Delta Shocks and Sin¬ 
gular Shocks, Nonlinear Theory of Generalized Eunctions, 
1999 

^ E. Yu. Panov and V. M. Shelkovich, S'-Shock waves as 
a new type of solutions to systems of conservation laws. 
Journal of Differential Equations, 2006 
R. Meppelink, S. B. Roller, J. M. Vogels, P. van 
der Straten, E. D. van Ooijen, N. R. Heckenberg, H. 
Rubinsztein-Dunlop, S. A. Haine, and M. J. Davis, Obser¬ 
vation of shock waves in a large Bose-Einstein condensate. 
Physical Review A, 2009 

R. J. Taylor, D. R. Baker, and H. Ikezi, Observation of 
collisionless electrostatic shocks. Physical Review Letters, 
1970 


Kurganov, A. and Noelle, S. and Petrova, G., Semidiscrete 
central-upwind schemes for hyperbolic conservation laws 
and Hamilton-Jacobi equations, SIAM J. Sci. Comput., 
2001 

J. P. Hirth and J. Lothe, Theory of Dislocations, 2nd ed., 
Krieger Publishing Company, 1982 

S. Limkumnerd and J.P. Sethna, Mesoscale theory of grains 
and cells: crystal plasticity and coarsening. Physical Re¬ 
view Letters, 2006 

A. Roy and A. Acharya, Finite element approximation of 
field dislocation mechanics, J. Mech. Phys. Solids, 2005 
D. A. Hughes, D. C. Chrzan, Q. Liu, and N. Hansen, Scal¬ 
ing of Mis orientation Angle Distributions, Physical Review 
Letters, 1998 

W. Choi, The Physics of Singular Dislocation Structures 
in Continuum Dislocation Dynamics, Thesis, Cornell Uni¬ 
versity, 2013 

Y. S. Chen, W. Choi, M. Bierbaum, J. P. Sethna, and S. 
Papanikolaou, Scaling theory of continuum dislocation dy¬ 
namics in three dimensions: Self-organized fractal pattern 
formation. International Journal of Plasticity, 2013 
S. Limkumnerd, Mesoscale Theory of Crains and Cells: 
Polycrystals & Plasticity Thesis, Cornell University 2007 
S. Limkumnerd and J. P. Sethna, Mesoscale Theory of 
Crains and Cells: Crystal Plasticity and Coarsening, Phys. 
Rev. Lett. 96, 095503 (2006). 

Yong S. Chen, Woosong Choi, S. Papanikolaou, and J. P. 
Sethna, Bending Crystals: Emergence of Fractal Disloca¬ 
tion Structures, Phys. Rev. Lett. 105, 105501 (2010). 
Woosong Choi, Yong S. Chen, Stefanos Papanikolaou, and 
James P. Sethna, Is dislocation flow turbulent in deformed 
crystals?. Computing in Science and Engineering, 14, 33 
( 2012 ). 



